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Abstract 

The thermodynamic driving force in the self-assembly of the secondary structure of a class 
of donor-acceptor oligorotaxanes is elucidated by means of molecular dynamics simulations of 
equilibrium isometric single-molecule force spectroscopy AFM experiments. The oligorotaxanes 
consist of cyclobis(paraquat-p-phenylene) rings threaded onto an oligomer of 1,5-dioxynaphthalenes 
linked by polyethers. The simulations are performed in a high dielectric medium using MM3 as 
the force field. The resulting force vs. extension isotherms show a mechanically unstable region in 
which the molecule unfolds and, for selected extensions, blinks in the force measurements between 
a high-force and a low-force regime. From the force vs. extension data the molecular potential of 
mean force is reconstructed using the weighted histogram analysis method and decomposed into 
energetic and entropic contributions. The simulations indicate that the folding of the oligorotaxanes 
is energetically favored but entropically penalized, with the energetic contributions overcoming the 
entropy penalty and effectively driving the self-assembly. In addition, an analogy between the 
single-molecule folding/unfolding events driven by the AFM tip and the thermodynamic theory of 
first-order phase transitions is discussed and general conditions, on the molecule and the cantilever, 
for the emergence of mechanical instabilities and blinks in the force measurements in equilibrium 
isometric pulling experiments are presented. In particular, it is shown that the mechanical stability 
properties observed during the extension are intimately related to the fluctuations in the force 
measurements. 

PACS numbers: 36.20.Ey, 82.37.Gk, 82.60.Qr 
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I. INTRODUCTION 



Unfolding and unbinding events in individual macromolecules can be computationally 
studied by means of "pulling" molecular dynamics (MD) simulations . These calcu- 

lations simulate recent single-molecule experiments in which a molecule (typically a protein 
or DNA) attached to a surface is mechanically unfolded by pulling it with optical tweezers 
or an atomic force microscope (AFM) tip attached to a cantilever 0, Is, [a 7|. Such ex- 
periments show stress maxima that have been linked to the breaking of hydrogen bonds, 
providing insight into the secondary and tertiary structure of the macromolecules, and the 
dynamical processes involved in muscle contraction, protein folding, transcription, shape 
memory materials and many others. 

In this paper we present simulations of single-molecule pulling for a class of donor-acceptor 
oligorotaxanes, as well as a detailed interpretation of the observed phenomenology. The pseu- 
dorotaxanes |8|, [Oj considered consist of a variable number of cyclobis(paraquat-]9-phenylene) 
tetracationic cyclophanes (CBPQT^"'") threaded onto a linear chain composed of three naph- 
thalene units linked by polyethers, with polyether caps at each end (see Fig. [T]). These 

are oligomeric analogues 

of polyrotaxanes previously developed by the Stoddart group pLl] . The relevant interactions 
that drive the folding of this class of molecules have, in the absence of applied stress, been 
recently characterized [12^. The interest in this paper is to investigate folding/unfolding 
pathways for these molecules under stress, and to determine, from a thermodynamic per- 
spective, what drives the folding. In this sense, the AFM pulling setup is a valuable tool 
since it allows inducing unfolding/refolding events at the single-molecule level while simul- 
taneously performing thermodynamic measurements. 

The AFM pulling of the rotaxanes is simulated using constant temperature MD in a 
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15| . We focus on the isometric 



high dielectric medium using MM3 as the force field 
version of these experiments in which the distance between the surface and the cantilever 
is controlled while the force exerted is allowed to vary. Closely related experimental efforts 
are currently underway in the Stoddart group. To bridge the several orders of magnitude 
gap between the pulling speeds that are experimentally employed (1 — 10^'^ /im/s) and 
those that are computationally accessible, in the simulations the unfolding events are driven 
slowly enough that quasistatic behavior is recovered and the results become independent of 
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the pulling speed. 

A central quantity of interest in this analysis is the molecular potential of mean force 
(PMF) 161] along the extension coordinate. The PMF is the Helmholtz free energy pro- 
file along a given reaction coordinate. As such, it succinctly captures the thermodynamic 
changes experienced by the molecule during the folding. Here, the PMF is reconstructed 
from the equilibrium force measurements using the weighted histogram analysis method 
(WHAM) 17, Q, [3]. Related approaches based on nonequilibrium force measurements 
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251 . l26l . |27[ . Further insight into the thermodynamic driving force behind the 
folding is obtained by decomposing the PMF into energetic and entropic contributions. 

In addition, an analogy between the phenomenology observed during the pulling of the 
rotaxanes and that expected for a system undergoing a first-order phase transition is pre- 
sented. In the process of exploring such an analogy, several general features of the isometric 
pulling experiments will be revealed. Specifically, the dependence of the force vs. exten- 
sion curves on the cantilever spring constant is characterized, and general conditions on the 
molecule and on the cantilever necessary for the emergence of mechanical instability and 
blinks in the force measurements are isolated. In particular, it is shown that the mechanical 
stability properties observed during the pulling are intimately related to the fluctuations 
in the force measurements. These results complement previous studies done by Kreuzer et 
al. 



in which the soft-spring and stiff-spring limit of the AFM pulling were identified with 
the Gibbs and Helmholtz ensembles for the isolated molecule, those of Kirmizialtin et al. |29 1 



in which the origin of the dynamical bistability observed in constant force experiments 
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and its relationship with the topography of the PMF was clarified, and those of Friddle 
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was analyzed for a simple model potential meant to represent bond-breaking (see Ref. 
for an analysis of the bond survival probability). The conditions isolated below apply to 
any system subject to equilibrium isometric pulling using AFM tips attached to harmonic 
cantilevers. 

This manuscript is organized as follows. In Sec. Ill Al the protocol employed to simulate 
the isometric pulling of the rotaxanes is described. Sections III Bl and III CI summarize, re- 
spectively, the strategy used to reconstruct the PMF from the force measurements and the 
procedure employed to decompose it into energetic and entropic contributions. Our main 
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results are presented in Sec. IIIII Specifically, in Sec. IIII Al the basic features observed dur- 
ing the mechanical unfolding of the oligorotaxanes are presented and the resulting PMF 
discussed. In Sec. IIIIBI the origin of the mechanical instabilities and blinks in the force 
measurements observed during the extension is clarified, and minimum conditions for the 
emergence of these effects are presented. Our main conclusions are summarized in Sec. IIV[ 

II. THEORETICAL METHODOLOGY 
A. Pulling simulations 

Unfolding of the rotaxanes is computationally studied by means of pulling MD simula- 
tions. The simulations are analogous to single-molecule force spectroscopy experiments in 
which a molecule attached to a surface is mechanically unfolded by pulling it with an AFM 
tip attached to a cantilever. The general setup of this class of experiments is shown in Fig. [2l 
The stretching computation begins by attaching one end of the molecule to a stiff isotropic 
harmonic potential that mimics the molecular attachment to the surface. Simultaneously, 
the opposite molecular end is connected to a dummy atom via a virtual harmonic spring. 
The position of the dummy atom is the simulation analogue of the cantilever position, and 
is controlled throughout. In turn, the varying deflection of the virtual harmonic spring mea- 
sures the force exerted during the pulling. The stretching is caused by moving the dummy 
atom away from the molecule at a constant speed. The pulling direction is deflned by the 
vector connecting the two terminal atoms of the complex. Since cantilever potentials are 
typically stiff in the direction perpendicular to the pulling, in the simulations the terminal 
atom that is being pulled is forced to move along the pulling direction by introducing ap- 
propriate additional harmonic restraining potentials. This type of simulations does not take 
into account any effects that may arise due to the interaction between the molecule and the 
AFM tip or the surface that cannot be accounted for by simple position restraints, nor any 
variations in the cantilever position due to thermal fluctuations or solvent-induced viscous 
drags. 

The potential energy function of the molecule plus restraints is of the form 



where UQ{r) is the molecular potential energy plus potential restraints not varied during the 
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simulation, r denotes the position of the atoms in the macromolecule, and 

VL[ar),t] = ^[ar)-mf (2) 

is the potential due to the cantilever. Here ^ (r) is the molecular end-to-end distance function, 
L = Lq + vI is the distance from the cantilever to the surface at time t, v is the pulling speed, 
and k the cantilever spring constant. The force exerted by the cantilever on the molecule at 
time t is given by 

F{t) = -VV, = ^ = -k lar) - m] , (3) 

where the gradient V is with respect to the [C, (r) — L] coordinate. In writing Eqs. ([2]) and ([3]) 
the vector nature of v(t), L(t), ^(r) and F(t) has been obviated since these quantities are 
collinear in the current setup. 

n 

The simulations are performed using Tinker 4.2 [32] for which a pulling routine was 
developed. The stretching is performed using a soft cantilever potential with spring constant 
k = 0.011 N/m = 1.1 pN/A. Since thermal fluctuations in the force depend on k as 6F ~ 
(where /3 = l/ksT is the inverse temperature), such a soft spring constant provides 
high resolution in the force measurements. The macromolecules are described using the MM3 
force field, which we have found adequately reproduces the complexation energies and the 



interactions responsible for the folding of the oligorotaxanes 



12| . As a simple model for the 



solvent, we use a continuum high dielectric medium in which the overall dielectric constant 
of the force field is set to that of water at room temperature (78.3). This is to be compared 

n 

with the MM3 dielectric constant for vacuum of 1.5 [33|. Any solvent effects that cannot be 
described using this continuum description are absent. The dynamics is propagated using 
a modified Beeman algorithm with a 1 fs integration time step, and the system is coupled 
to a heat bath at 300 K using a Nose-Hoover chain as the thermostat. Initial -minimum 
energy- structures for the pulling simulations were obtained using simulated annealing, as 
described in Ref. [l^- These initial structures were allowed to equilibrate thermally for 1 ns 
and subsequently stretched. The simulations do not include counterions for the tetracationic 
cyclophanes since, at 300K and for the high dielectric medium employed, their interaction 
with the main molecular backbone is small and their effect on the folding negligible. 
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B. Reconstructing the PMF using WHAM 



Under reversible conditions, knowledge of the force exerted during the pulling immedi- 
ately yields the associated change in the Helmholtz free energy A for the molecule plus 
cantilever. This is because the change in A is determined by the reversible work exerted 
during the pulling 

AA^£(^).d.'^-il„|W. (4) 

where (F) l is the average force at extension L, and 

Z{L) = J drexp[-/3f/z.(r)] (5) 

is the configurational partition function of the molecule plus cantilever. Equation (jl]) follows 
from standard thermodynamic integration considerations [3], and assumes that at each 
point during the extension the state of the molecule plus cantilever potential is well described 
by a canonical ensemble. This property is expected for a system in weak contact with a 
thermal bath and is enforced in the simulations by the nature of the thermostat employed. 

From the force versus extension measurements it is also possible to extract the molecular 
potential of mean force (PMF) 0(^) as a function of the end-to-end distance ^ by properly 
removing the bias due to the cantilever potential. The quantity 0(^) is the Helmholtz 
free energy profile along the coordinate ^ for the isolated (cantilever-free) molecule, and is of 
central interest since it succinctly characterizes the thermodynamics of the unfolding. Below 
we summarize how to estimate this quantity from the force measurements employing the 
weighted histogram analysis method (WHAM). 



The PMF is defined by |34| 



m = 0(r ) - ^ In 

where ^* and 0(^*) are arbitrary constants and 
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. /d..[^-e(.)]expM^o(.)] ^ M) . ,Si, - (7) 
j dr exp[-/3f/o(r)J Zq 



Here, 



Zo = Jdr exp[-(3Uo{r)] (8) 
is the configurational partition function of the molecule, and 

ZoiO = J dr5[e-e(r)]exp[-/3f/o(r)]. (9) 



In the context of the pulhng experiments PoiO is the probabihty density that the molec- 
ular end-to-end distance function ^(r) adopts the value ^ in the unbiased (cantilever-free) 
ensemble. The quantity Po(0 determines the PMF up to a constant and can be estimated 
from the force measurements, as we now describe. 

At this point it is convenient to discretize the time variable during the pulling experiment 
into M steps: ti, . . . ,ti, . . . ,tM- The potential of the system plus cantilever at time U is 
Ui = Uo + Vi, where Vi = |[,^(r,tj) — L(tj)]^ is the bias due to the cantilever potential. In the 
2-th biased measurement knowledge of the force Fi and L{ti) gives the molecular end-to-end 
distance. The probability density of observing the value ^ at this time is given by 

j dr exp[-/3Ui{r)\ Zi 
where Zi = J dr exp[— j3Ui{r)] is the configurational partition function for the system plus 
cantilever at the i-th extension, and 

Z,(0 = J dr5[^ - e(r)] exp[-/5f/,(r)]. (11) 

Knowledge of Pi{^) allows the unbiased probability density Po(0 to be reconstructed since, 
by virtue of Eqs. and ffTU]) . these two quantities are related by 

Po(0=exp[W.(0]|^P.(0, (12) 

where we have exploited the fact that Zq{^) = exp[+/9Vi(^)]Zj(^). In the pulling experiments, 
for each L{ti) only a few force measurements are typically performed and the probability den- 
sity Pi{^) is not properly sampled. Nevertheless, the experiments do provide measurements 
at a wealth of values of L{ti) that can be combined to estimate Po(0- 

M „ 
PoiO = y2'^^eXp[+mO]l^P^{0, (13) 

where the Wi are some normalized (^^ Wi = 1) set of weights. In the limit of perfect 
sampling any set of weights should yield the same po(0- practice, for finite sampling it is 
convenient to employ a set that minimizes the variance in the po(0 estimate from the series 
of independent estimates of biased distributions. Such a set of weights is precisely provided 
by the WHAM prescription 0, 3, 19, 35 1 



where Oi = 1 + 2rj is the statistical inefficiency 35(] and the integrated autocorrelation 



time 
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MO = ^r.:^L^ (15) 



351]. The difficulty in estimating the Tj for each measurement generally leads to 
further supposing that the QiS are approximately constant and factor out of Eq. (fill) , so 
that 

Neglecting the gis from Eq. (ITS!) does not imply that the resulting estimate of Po(0 
incorrect; but simply that the weights selected do not precisely minimize the variance in the 
estimate. Experience with this method indicates that if the QiS do not differ by more than 
an order of magnitude their effect on 0(^) is small [l^. 

The computation of the PMF then proceeds as follows. Suppose that Ni force measure- 
ments are done for each extension i. Knowledge of the force {j = 1, . . . ,Ni) and of L{ti) 
gives the molecular end-to-end distance in each of these measurements. The numerator 
in Eq. (fT5|) is then estimated by constructing a histogram with all the available data, 

i=l i=l j=l * ^ 

where A,^ is the bin size, and C| (^) = 1 if ^ ['C ~ + ^^/2) and zero otherwise. 

Estimating the denominator in Eq. (fT5|) requires knowledge of the Zj's, the configurational 
partition functions of the system plus cantilever at all extensions considered. There 

are two ways to obtain these quantities. The most direct one is to employ the Helmholtz 
free energies for the system plus cantilever obtained through the thermodynamic integration 
in Eq. (jl]), as they determine the Zj's up to a constant multiplicative factor. Alternatively, 
it is also possible to determine the ratio of the configurational partition functions between 
the i-th biased system and its unbiased counterpart using Po{C,)- 

Y^=J^^ ^M-mOMO- (17) 

Equations (ITS!) and (ITTI) can be solved iteratively. Starting from a guess for the Zi/Zo, Po{^) 
is estimated using Eq. f|T5|) and normalized. The resulting po(0 is then used to obtain a 
new set of Zi/Z^ through Eq. (ITTI) . and the process is repeated until self-consistency. This 
latter approach is the usual procedure to solve the WHAM equations. Note, however, that 
this procedure is not required if reversible force vs. extension data is available. 

The reconstruction of the PMF from the force vs. extension measurements using WHAM 
is simple to implement. Further, the procedure is versatile in the sense that it makes no 
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assumption about the stiffness of the cantilever potential and permits the combination of 
data from several pulling runs. In addition, as we describe below, it can also be employed 
to decompose the PMF into entropic and energetic contributions from data readily available 
from an MD run. 



C. Energy-entropy decomposition of the PMF 

Insight into the thermodynamic driving force responsible for the folding of oligorotaxanes 
can be obtained by decomposing the changes in the PMF into entropic S{C,) and potential 
energy t/o(0 contributions, 

m - = [Uoio - uoin] - t m) - s{o] ■ m 

Here, 

m) = ^/ dr(5[e-e(r')]t/o(r)exp[-/5[/o(r)] (19) 
is the average potential energy when the end-to-end distance adopts the value ^, and 

TS{i) = j 5[i - i{r)]p\n p (20) 

is the corresponding configurational entropy, where p = exp[— /3f/o(r)] and VIq is a 
constant irrelevant factor with the same dimensions as Zq[E,)- Equation (|T8|) follows from 
the above definitions. 

The energy-entropy decomposition of the PMF can be obtained by estimating Uo{S,) di- 
rectly and employing the previously reconstructed 0(^) to obtain TS{^) using Eq. (fT5]l . Here, 
Uq{S,) is estimated by combining all measurements performed during the pulling simulations 
using WHAM. The measurable quantity in the simulation in the presence of the bias due to 
the cantilever is: 

mrm-iir)^),^ 5^^m-iirWM-mAr)] ^ (21) 

the average molecular potential energy when the end-to-end distance adopts the value ^ 
and the cantilever extension is L{ti). This average in the biased ensemble is related to the 
unbiased average f/o(0 by: 

UoiO = ^M+mO]^{Uoir)5[^ - ar)]h. (22) 



As in the estimation of po{^), all data collected during pulling can be combined with appro- 
priate weights Wi to estimate f/o(0' 

M 

m) = J2''^^M+mO]^AUo{r)6[^-ar)]h. (23; 
If one adopts the weights Wi that define the WHAM prescription [Eq. (HM ]. then 



where we have assumed identical statistical inefficiencies for all i. Last, recalling the defini- 
tion of Po(0 = Zo{^)/Zo [Eq. (171)] and its WHAM prescription [Eq. ( ITSi) ]. one arrives at a 
useful expression for Uo{^), 

Equation (125|) provides a mean to calculate the potential energy as a function of the end- 
to-end distance by combining all the data obtained during the pulling. The remarkable 
simplicity of Eq. (125|) is noteworthy: in order to obtain t/o(0 just needs to generate 
histograms for the the values of ^ observed using all data collected during the pulling and 
then, within each bin, perform a simple average of the internal energy of the configurations 
that fall into it. It does not require knowledge of the properties of the bias potential. 

This method of performing an energy-entropy decomposition of the PMF requires knowl- 
edge of ?7o(0 the configurations observed during pulling. This quantity, although readily 
available in an MD run, is not experimentally accessible. As a result, in experiments an 
energy-entropy decomposition of the PMF requires estimating the changes in entropy by 
performing the pulling at varying temperatures. 

III. RESULTS AND DISCUSSION 

A. Phenomenology of the pulling experiments 

1. The approach to equilibrium 

One of the challenges in simulating AFM pulling experiments using MD is to bridge the 
large disparity between the pulling speeds v that are employed experimentally and those 
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that can be accessed computationally. While typical experiments often use f ~ 1 — 10~^ 
/im/s, current computational capabilities require pulling speeds that are several orders of 
magnitude faster. Here, this difficulty is circumvented by striving for pulling speeds that are 
slow enough that reversible behavior is recovered. Under such conditions the results become 
independent of f , and the simulations comparable to experimental findings. 

Consider the pulling of the [3]rotaxane shown in Fig. [1] (the simulation details are spec- 
ified in Sec. Ill A|) . Figure |3] shows the effect of decreasing the pulling speed on the force 
vs. extension characteristics. In the simulations, the system is first extended (black lines) 
to a given L and then contracted (gray lines). For pulling speeds (and equal retracting 
speeds) between 1 — 0.01 A/ps (panels A to C) hysteresis and other nonequilibrium effects 
during the pulling are unavoidable. However, for v = 0.001 A/ps = 10^ /xm/s the system 
behaves reversibly, and the F-L curves obtained during extension and contraction essen- 
tially coincide. Note the second law of thermodynamics at play in the simulations: the work 
required to stretch the molecule -the area under the curve- decreases with the pulling speed 
until reversible behavior is attained. We have observed that pulling speeds of ~ 10"'^ A/ps 
are generally sufficient to recover reversible behavior for the unfolding in a high dielectric 
medium of the [n]rotaxanes considered. For simulations in vacuum slower pulling speeds are 
required. 

It is worth noting that this seemingly simple phenomenological behavior is not recovered 
when a Berendsen thermostat is employed instead of a Nose- Hoover chain. In fact, we have 
observed that the Berendsen thermostat leads to a spurious violation of the second law 
during the pulling, presumably because the ensemble generated by it does not satisfy the 
equipartition theorem 36 1. 



2. Reversible unfolding, mechanical instability and blinks in the force measurements 

Figure H] details the force exerted (upper panel) and the molecular end-to-end distance 
(lower panel) during the reversible pulling and contraction of the [3]rotaxane. Snapshots 
of typical structures encountered during the pulling are included in Fig. [51 As the system 
is stretched the oligorotaxane undergoes a conformational transition from a folded globular 
state to an extended coil. In the process, the force initially increases approximately linearly 
with L, then drops, and subsequently increases again. The drop in the force is due to the 
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unfolding of the rotaxane during the stretching. 

Figure [6] shows the change in the Helmholtz free energy A for the molecule plus cantilever 
obtained from integrating directly the F-L curve in Fig. m The net change in free energy for 
the complete thermodynamic cycle is zero, as expected for a quasistatic process. The free 
energy profile for the system plus cantilever has three distinct regions, labeled I-III in the fig- 
ure. Regions I and III correspond to the folded and extended state, respectively, while region 
II is where the folding/unfolding event occurs. From a thermodynamics perspective 37|, 
regions I and III are mechanically stable phases since A is a convex function of the exten- 



sion L, i.e. (^firj = y dh ) — contrast, region II, where the unfolding occurs, is 

mechanically unstable since A is a concave function of L and hence ^f|4^ = { ^^ql^ j < 0- 
The dynamical behavior of the system in the unstable region is illustrated in Fig. [71 which 
shows the evolution of the radius of gyration Rg and the molecular end-to-end distance when 
L = 70.0 A. In this region, the [3]rotaxane undergoes transitions between a folded globular 
state, partially folded structures, and an extended coil. The right panels in Fig. [7] show the 
probability density of the distribution of and ^ values obtained from a 20 ns trajectory. 
The system exhibits a clear dynamical bistability along the ^ coordinate, and at least a 
tristability along R^. Figure [HI shows some representative structures encountered during this 
dynamics. The bistability along the end-to-end distance leads to a blinking in the force 
measurements from a high force to a low force regime during the pulling (cf. the unstable 
region in Fig. H]). The blinking in the force measurements in the isometric experiments 
is analogous to the blinking in the molecular extension observed during the constant force 
reversible pulling of single RNA molecules ^] . Note that for the cantilever stiffness employed 
this multistable behavior in the radius of gyration, or the end-to-end distance, is not observed 
when L is set to be in the stable regions (I or III) of the free energy profile. The origin of 
the bistability along the end-to-end distance, the multistability in the radius of gyration, 
and its relation to the thermodynamic instability in the pulling is discussed in Sec. IIIIBI 



3. The PMF and the thermodynamic driving force in the folding 

Figure [9] shows the potential of mean force along the end-to-end distance C, for the [3]ro- 
taxane reconstructed from the force measurements as described in Sec. Ill B| as well as its 
decomposition into entropic and potential energy contributions. The thermodynamic native 
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state of the oligorotaxane, i.e. the minimum in 0, corresponds to a globular folded structure 
with ^0 = 7 A. The structure of the molecular free energy profile as a function of ^ is similar 
to that observed for the molecule plus cantilever as a function of L. Namely, it consists of a 
folded and an unfolded stable branch where </)(^) is a convex function, and a concave region 
around ^ = 44 — 52 A where the conformational transition occurs. As discussed in Sec. IIII Bt 
this characteristic structure of the PMF is responsible for many of the interesting features 
observed during pulling. Figure [10] shows the probability distribution for i?^ at values of 
^ in the three main regions of the PMF, each collected during a 20 ns trajectory. In the 
stable regions = 7.0 and 64.9 A ), the evolution of the radius of gyration reveals only one 
stable structural state. By contrast, in the unstable region = 47 A) the distribution of 
R^g values indicates that there is a inherent bistability in the molecular potential along the 
R^g coordinate. Such molecular bistability gives rise to the region of concavity in the PMF. 

What drives the folding of the oligorotaxanes from a thermodynamic perspective? The 
energy-entropy decomposition of the PMF indicates that the folding of the rotaxane is 
energetically favored but entropically penalized, with the energetic contributions overcoming 
the entropy penalty and effectively driving the self-assembly. Further, while the potential 
energy and the entropy of the molecule show considerable changes during the unfolding, these 
two effects largely cancel one another, leading to only modest changes in the Helmholtz free 
energy. Note that while is a monotonic function of ^ for ^ > ^q, the energy and entropy 
contributions show a clear nonmonotonic dependence on ^. For instance, there is a reduction 
in the system's entropy for ^ both just before the conformational transition and for large 
extensions due to a reduction in the available degenerate conformational space. 

It is important to stress that the PMF and its decomposition into energy and entropy 
contributions correspond to a molecule in which the terminal atoms are constrained to move 
along the pulling coordinate. As a consequence, any additional contributions that may arise 
by relaxing these constraints are simply not manifest in this picture of the folding. 

4- Dependence of the PMF on the number of threaded rings 

Figure [TT] compares the PMF extracted from force measurements for the [3]rotaxane, 
and the associated [2]rotaxane and bare molecular thread. In the [2]rotaxane system the 
CBPQT^+ring encircles the central naphthalene unit of the underlying chain. The PMF 
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in all cases is qualitatively similar: it consists of two stable convex regions, and a unstable 
concave region where the unfolding occurs. Further, for all species the folding along is 
energetically driven and entropically penalized (not shown). Note that, in the high dielectric 
medium employed, the number of threaded rings does not have a dramatic effect on the net 
change in the free energy undergone during folding nor, as expected, on the elastic properties 
of the extended state. It does, however, have an appreciable effect on the stability and 
elasticity of the folded conformation. Specifically, as the number of rings is increased the 
range of ^ values for which the folded conformation is stable increases. That is, a larger L 
is required during the pulling to unfold the rotaxanes with respect to the bare chain. 



B. Interpretation of the pulling experiments 

1. Analogy with first-order phase transitions 

At first glance, the phenomenological behavior observed by the molecule plus cantilever 
during the extension is analogous to that expected for a system undergoing a first-order 



phase transition [37|]. Namely, during the extension the combined system goes from one 
stable thermodynamic phase to another (Regions I and III in Fig. [6]) by varying an externally 
controllable parameter, in this case L. In the transition region the Helmholtz free energy 
of the combined system is concave and, hence, fails to satisfy the thermodynamic stability 
criteria (Region II in Fig. [6]). As a consequence, the observed F-L isotherm (Fig. H]) exhibits 
a region of mechanical instability where ^^^^j < 0. These observations largely parallel 
the behavior of unstable P-V isotherms for a van der Waals fiuid. 

In addition, the phenomenology seems to indicate that close to the transition region the 
Helmholtz potential for the molecule plus cantilever is bistable along the end-to-end distance, 
and has the characteristic form exhibited during a first-order phase transition schematically 
shown in Fig. [121 For L G Region I the stable global minimum corresponds to the molecule 
in a folded state, and as L is increased the equilibrium state shifts from one local minimum to 
the other. For L G Region II the two minima become approximately equal and a transition 
from the folded state to the unfolded state occurs. Beyond this extension the unfolded 
molecular phase becomes the global potential minimum and the extended phase becomes 
absolutely stable. This interpretation is consistent with the observed behavior in which the 
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system blinks between a folded and an extended state around the transition region, recall 
Fig. [3 

Note, however, that this analogy can never be completely faithful since the modeled 
single-molecule pulling experiment does not deal with a macroscopic system. This contrasts 



with the elastic behavior observed for polymer chains |38l . |39| , and leads to salient differences 
between the two processes. First, while in a true first-order phase transition the physical 
system exhibits coexistence between the stable phases, in this single-molecule version no 
coexistence is possible. At the transition point, and for the value of k employed, the system 
instead exhibits frequent jumps between the stable phases (cf. Fig. [7j). Coexistence is re- 
placed by ergodicity in this single-molecule manifestation of bistability. Second, the variable 
that is varied during the extension L is not a true thermodynamic variable since it is neither 
extensive nor intensive. Third, the system considered has no clear thermodynamic limit 
in which a true discontinuity of the derivatives of the partition function can develop and, 
further, the system is far from the thermodynamic limit since the fluctuations observed in 
the F-L curves are c omp arable to the average value, leading to a nonequivalence between 



statistical ensembles 
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411]. Additionally, the mechanical properties measured during the 
pulling are for the molecule plus cantilever, and hence the extension behavior is expected to 
depend on the nature of the cantilever. 

Nevertheless, the similarities between the two processes are still intriguing and, in view 
of these observations, it is natural to ask: (i) what is the origin of the mechanical instability 
during the puUing?; (ii) how does the dynamical bistability along the end-to-end distance 
arise?; (iii) how does the nature of the cantilever affect the observed phenomenological 
behavior; and (iv) which features of the observed behavior arise due to the molecule, and 
which due to the cantilever? These questions are addressed in the following sections. 

The key quantity in the remainder of this analysis is the PMF along the end-to-end 
distance [Eq. ([6])]. Its utility relies on the fact that the configurational partition function of 
the molecule plus cantilever at extension L [Eq. ([5])] can be expressed as 

Z(L) = Jd^ exp{-mO + VLiO]}. (26) 

where we have neglected the constant and L-independent multiplicative factor that arises 
in the transformation since it is irrelevant for the present purposes. The above relation 
implies that the extension process can be viewed as thermal motion along a one-dimensional 
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effective potential determined by the PMF 0(^) and the bias due to the cantilever Vl(^): 



ULiO = m + VLiO- (27) 

Further, by means of the PMF it is possible to estimate the force vs. extension characteristics 
for the composite system for any value of the force constant k. This is because the average 
force exerted on the system at extension L can be expressed as: 

l d\n[Z{L)/Zo] _ /deqfflexp{-/3[0(O + r,(O]} 
^ f3 dL Jd^exp{-/3[0(O + l^L(O]} ' ^ ' 

where, for convenience, we have introduced Zq in the logarithm just to make the argument 
dimensionless. Since </)(^) is a property of the isolated molecule it is independent of k. Hence, 
Eq. (!28|) can be employed to estimate the F-L isotherms for arbitrary k. 



2. Dependence of the F-L isotherms on the cantilever spring constant 

Figure [13] shows representative F-L curves for the [3]rotaxane estimated using Eq. (|28|) 
for different cantilever spring constants, as well as the ratio between the thermal fluctuations 
and the average in the force measurements. Note that in all cases the system is far from 
the thermodynamic limit since the force fluctuations are comparable in magnitude to the 
average. Further, the fluctuations increase with the spring constant and, for example, in the 
case of = 5^0 (where k^ = 1.1 pN/A) the region of instability in the F-L curve would be 
mostly masked by the fluctuations. 

Note that the instability properties of the F-L isotherms depend intricately on the can- 
tilever spring constant employed. This dependence is succinctly conveyed in Fig. [H] which 
shows the values of the extension and the average force at the critical points in the F-L 
isotherms. In the figure, F~^ and L"*" (or F~ and L~) denote the values of the force and the 
cantilever extension when the F-L curve exhibits a maximum (or minimum). These values 
enclose the unstable region in the F-L isotherm. In turn, is the average end-to-end dis- 
tance at the critical points which is approximately independent of k except for very small k. 
The critical forces F~^ and F~ show a smooth and strong dependence on the cantilever spring 
constant. For large k the system persistently shows a region of instability in the isotherms. 
However, as the cantilever spring is made softer, F~^ and F^ approach each other and for 
small k the mechanical instability in the F-L isotherms is no longer present. 
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3. Disappearance of the mechanical instability in the soft spring limit 



The disappearance of the instabihty regions in the F-L measurements when the system 
is pulled with a soft cantilever is a general feature of the extension of any molecular system 
provided that no covalent bond breaking occurs during the process. To see this, consider 
the soft-spring approximation of the configurational partition function of the system plus 
cantilever [Eq. fl26l) ]: 

^ exp(-^fcLV2) r d^exp[-/30(O]exp(/?fcLO = exp(-/3A;LV2)(exp(/?fcLO), (29) 

where the notation (/) stands for the unbiased (cantilever-free) average of /. In writing 
Eq. (129!) we have supposed that in the region of relevant ^ (in which the integrand is non- 
negligible) k^"^ /Vl{^) <^ 1 and, hence, that the cantilever potential is well approximated 
by Vl{C} ~ — kL^. This approximation is valid provided that that no bond breaking 
is induced during pulling and permits the introduction of a cumulant expansion 42] in the 
configurational partition function. Specifically, the average {exp{PkL^)) can be expressed 
as: 

ieMPkLO) = f: = exp jf; ^-n(0| (30) 

n=0 ■ Ln=l ■ j 

where Kn(0 is the n-th order cumulant. In view of Eqs. (1281) . (!29|) and (!30!) the slope of the 
F-L curves can be expressed as: 

^-*-^*=E^-«(o. (31) 

n=0 

It then follows that to lowest order in k, 

^^k>0. (32) 

That is, no unstable region in the F-L curve can arise in the soft-spring limit, irrespective of 
the specific form of 0(0- Instabilities, however, can arise in higher orders in the expansion. 
For instance, to second order in k instabilities arise provided that the cantilever spring 
constant satisfies: 



where K2{0 = (e) - iO'- 
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4- Emergence of the mechanical instability in the stiff spring limit 



A mechanically unstable region in the F-L curves can only develop if the PMF of the 
unbiased molecule by itself has a region of concavity. To understand how this result arises 
consider the configurational partition function of the molecule plus cantilever at extension 
L [Eq. (!26|) ] for large k. In this regime, most contributions to the integral will come from 
the region where ^ = L. Consequently, exp[— can be expanded around this point to 



give: 



exp[-/30(O]=exp[-/50(L)] 



i-(3A,{m-L)-^A,{m-Lr+.. 



(34) 



where 



A,{L) 



dm 

dL 



...,^^-.(^)^ 



(35) 



Here, the notation is such that = 



. Introducing Eq. into Eq. (!26ll . 



integrating explicitly the different terms, and performing an expansion around 1/ A; = one 
obtains: 



Z{L) 



27r 




kj3 



exp[-/?<^(L)] 



1-^^(L) + 0(1/P) 



(36) 



Using Eq. (!28ll . the derivative of the force can be obtained from this approximation to the 
partition function. To lowest order in 1/ A; it is given by: 



dL dL^ 213k dL"^ 
from which it follows that, to zeroth order in 1/k, 

d{F)L 



+ Oil/k). 



(37) 



(38) 



dL dL^ 

That is, an unstable region in the force vs. extension requires a region of concavity in the 
PMF. 



5. Relation between force fluctuations and mechanical stability 

The mechanical stability properties observed during the AFM pulling experiments are 
intimately related to the fluctuations in the force measurements. To see this, note that the 
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average slope of the F-L curves can be expressed in terms of the fluctuations in the force: 



dL \ dL 



- (3[{F'), - {F)i] = k - f3[{F'), - {F)i], (39) 

L 

where we have employed Eq. fl28l) . The sign of determines the mechanical stabil- 

ity during the extension and hence Eq. (!39|) relates the thermal fluctuations in the force 



measurements with the stability properties of the F-L curves. Specifically, for the stable 

dL 



branches for which ^^D^ > the force fluctuations satisfy 



{F')l - {F)l < ^. (40a) 

In turn, at the critical points where the derivative changes sign the force fluctuations satisfy 
a strong constraint: 

{F')l - {F)l = ^. (40b) 

Last, in the unstable branches the force fluctuations are larger than in the stable branches 
and satisfy the inequality 

{F')l - {F)l > |. (40c) 

These inequalities can be employed to determine critical points in the force by measuring the 
fluctuations even in situations where the fluctuations mask the presence of critical points. 

Figure [T5] illustrates these general observations in the specific case of the pulling of [3]ro- 
taxane. As shown, for k = kQ and k = 2kQ there is a region where the force fluctuations 
become larger than k/(3 and, consequently, unstable behavior in F-L develops. For k = 0.4/co 
or less the fluctuations in the force are never large enough to satisfy Eq. fl40cl) and no critical 
points in the F-L develop, as can be confirmed in Fig. [TH 



6. Origin of the blinking in the force measurements 

The observation of force measurements that blink between a high-force and a low force 
regime (recall Fig. [7j) requires the composite system to be bistable along the end-to-end 
distance for some L, i.e. the effective potential Ul{0 [Eq- (1271) ] must have a double minimum. 
Since the PMF of the molecule by itself is not bistable along C, (see Fig. [9]), the bistability 
must be introduced by the cantilever potential. Figure [16] shows Ul{0 selected L for 
the [3]rotaxane. For L in the stable regions of the free energy (L = 40 A and L = 90 
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A) the effective potential exhibits a single minimum along the ^ coordinate. However, for 
L = 70 A a bistability in the potential develops. The secondary minimum is the cause for 
the blinking in the force measurements observed at this L (cf. Fig. [7]). Further, the relative 
energy between the two induced minima around the transition region can be manipulated 
by varying L. 

What are the minimum requirements for the emergence of bistability along ^? A necessary 
condition is that the effective potential Ul{C,) is concave for some region along ^, that is: 

^^ = ^ + ^<0- (41) 

For Eq. fj4T]) to be satisfied it is required that both: (i) the PMF of the isolated molecule 
has a region of concavity where < 0, and (ii) the cantilever employed is sufficiently 

soft such that 



k<-^A-^]. (42) 



for some ^. Equation (l42l) imposes an upper bound on > for bistability to be observable. 
If k is stiff the inequality would be violated for all ^ and bistability would not be manifest. 
Note, however, that there is no lower bound for k that prevents bistability along ^. This 
is in stark contrast with the behavior of the mechanical instabilities which disappear in the 
soft-spring limit and persist when stiff springs are employed (see Sees. 1111 B 211111 B 51) . 

Figure [17] illustrates these general observations during the pulling of the [3]rotaxane. 
Specifically, the figure shows, for cantilevers of varying stiffness, the probability density of 
observing the value F in the force measurements Pl{,F) (upper panels) and the associated 
probability density Pl{C) that the molecule adopts the extension ^ (lower panels) when the 
cantilever is fixed at L. For soft springs the system satisfies inequality (1421) for a range of 
L values and there is a clear bistability in both probability distributions. As k is increased 
this bistability becomes less prominent and for k = 5/co is no longer observed. 

7. Final remarks 

What then should be the picture of the pulling? Clearly, a picture like the one shown in 
Fig. [121 although appropriate to describe the phenomenology observed in Sec. IIII Al around 



the transition region (and the basis of an influential model of molecular extensibility 43|). 



is not accurate since the molecule itself does not show a double well structure along ^ in 
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the free energy. Note, however, that the molecular free energy is inherently bistable along 
some natural unfolding pathway. If one were able to determine such a pathway, a picture 
reminiscent to the one shown in Fig. [T2] would emerge. However, such natural unfolding 
pathway is not necessarily accessible during the pulling experiment. 

The analogy between first-order phase transition and the isometric single-molecule pulling 
experiments, even when suggestive, is limited by the fact that the qualitative features ob- 
served during pulling strongly depend on the cantilever constant employed. This is because 



the properties that are measurec 
However, as previously noted 



are those of the combined cantilever plus molecule system. 
, in the stiff spring limit the measured F-L curves resemble 
those which would be obtained from a calculation in the Helmholtz ensemble of the isolated 
molecule. In this limit, the instability behavior in the force vs. extension will only be due to 
molecular properties [recall Eq. fl38|l ]. and the associated bistability is that observed along 
the radius of gyration when ^ is fixed in the concave region of the PMF (Fig. [TOl) . It is 
in this limit that the phenomenology observed during the isometric single-molecule pulling 
is closest to a first-order phase transition. Experimentally, this limit is impractical since 
employing a stiff spring reduces the sensitivity in the force measurements. 



IV. CONCLUSIONS 



In this contribution we have simulated the equilibrium isometric AFM pulling of a series 
of donor-acceptor oligorotaxanes by means of constant temperature molecular dynamics in 
a high dielectric medium using MM3 as the force field. For this system, hysteresis and 
other nonequilibrium effects that may arise during pulling can be overcome by reducing the 
pulling speed. The resulting equilibrium force vs. extension isotherms show a mechanically 
unstable region in which the molecule unfolds and, for selected extensions, blinks in the 
force measurements between a high-force and a low-force regime. 

From the force vs. extension data, the Helmholtz free energy profile for the oligorotaxanes 
along the extension coordinate was reconstructed using the weighted histogram analysis 
method, and decomposed into energetic and entropic contributions. The simulations reveal 
an unfolding pathway for the oligorotaxanes, and indicate that the folding is energetically 
driven but entropically penalized. Even when the energy and entropy contributions to the 
potential of mean force can vary widely along the extension, these effects largely cancel one 
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another leading to only modest changes in the free energy profile. Further, we have observed 
that, for the rotaxanes studied, increasing the number of threaded rings stabilizes the folded 
conformation, making it resilient to unfolding over a wider range of end-to-end distances. 
As expected, the elastic properties of the extended conformation are not affected by the 
number of threaded rings. 

The potential of mean force of the rotaxanes exhibits only one minimum along the end-to- 
end distance. It consists of two stable convex regions characterizing the folded and unfolded 
conformations and a region of concavity where the unfolding occurs. This thermodynami- 
cally unstable region arises because the molecule is inherently bistable. The inherent molec- 
ular bistability is not resolved along the end-to-end coordinate (cf. Fig. ^ but, at least in 
the cases considered, it is manifest along the radius of gyration (recall Fig. [TU]) . 

Clearly, modifying the nature of the solvent employed during the pulling can have a 
quantitative effect on the determined free energy profile since the solvent can vary the 
relative stability of the different molecular conformations encountered along the extension. 
Nevertheless, except for very good solvents, the qualitative features of the topography of the 
potential of mean force described above are expected to remain intact. Also, hydrodynamic 
effects that may appear when considering an explicit solvent (e.g. drag forces acting on the 
AFM cantilever due to viscous friction with the surrounding solvent [44]) are immaterial for 
the equilibrium pulling. 

In addition, general conditions on the molecule and on the cantilever stiffness required for 
the emergence of mechanical instability and blinking in the force measurements in isometric 
pulling experiments have been presented. Specifically, we have shown that for such effects 
to arise the potential of mean force of the molecule needs to exhibit a region of concavity 
along the end-to-end distance. Further, the stiffness of the cantilever employed needs to 
satisfy specific constraints. For the blinks in the force measurements to arise the cantilever 
spring constant needs to be sufficiently soft so that Eq. ( |42|) is satisfied. In turn, for the 
mechanical instability to be observable the fluctuations in the force measurements have to 
be sufficiently large to satisfy Eq. f l40cp for some L. In practice, this implies that in the 
stiff-spring limit the mechanical instability can emerge, while in the soft spring limit the 
mechanical instability decays regardless of the details of the molecular system, provided 
that no covalent bond-breaking occurs during pulling. 

In this light, it becomes apparent that the analogy presented between first-order phase 
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transitions and the mechanical properties observed during the extension is hmited by the 
fact that the observed behavior depends on the cantilever constant employed. The analogy 
is the closest in the stiff-spring limit where the mechanical properties recorded are those 
that would have been obtained in the Helmholtz ensemble of the isolated molecule. 
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FIG. 1: Structure of the [3]rotaxane. Oxygen atoms are depicted in red, the naphthalene units 
in black, and the polyether carbons in orange. The CBPQT'^'^rings are depicted in blue with the 
pyridinium N"*" ions in brown. 



26 




FIG. 2: Schematic of an isometric single-molecule force spectroscopy experiment using an atomic 
force microscope (AFM). In it, one end of the molecular system is attached to a surface and the 
other end to an AFM tip attached to a cantilever. During the pulling, the distance between the 
surface and the cantilever L{t) = Lq + vt is controlled and varied at a constant speed v. The 
deflection of the cantilever from its equilibrium position measures the instantaneous applied force 
on the molecule by the cantilever F{t) = —k[S^{t) — L{t)], where k is the cantilever spring constant 
and (^(t) the fluctuating molecular end-to-end extension. 
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FIG. 3: The approach to equilibrium in the pulling simulations. The figure shows force vs. exten- 
sion profiles for the [3]rotaxane immersed in a high dielectric medium at 300K for different pulling 
speeds v of (A) 10°, (B) 10~^, (C) 10""^ and (D) 10~^ A/ps. The harmonic cantilever employed has 
a soft spring constant oi k = 0.011 N/m, L is the distance between the surface and the cantilever 
and F the instantaneous applied force. The black lines correspond to the extension process; the 
gray ones to the contraction. 
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FIG. 4: Time dependence of the force F and the molecular end-to-end distance ^ dm'ing the 
pulling of the [3]rotaxane under equilibrium conditions {v = 0.001 A/ps, k = 0.011 N/m). Typical 
structures (labels 1-7) observed during the extension are shown in Fig. [5j 
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FIG. 5: Snapshots of the [3]rotaxane during its extension. The numerical labels shown here are 
employed in Figs. I4|6I and [H 
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FIG. 6: Changes in the Helmholtz free energy A of the [3]rotaxane plus cantilever obtained from the 
data shown in Fig. [H The blue line corresponds to the pulling and the red one to the contraction. 
The dotted lines provide an estimate of errors in the thermodynamic integration due to force 
fluctuations at each pulling step. The difference in the degree of convexity of regions I and III is 
evidenced by extrapolating the data in each region outside of its domain through fitting to quartic 
polynomials (black dashed lines). The labels correspond to the structures shown in Fig. [5l 
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FIG. 7: Time dependence and probability density distribution of the radius of gyration (Rg) and 
the end-to-end molecular extension ^ when the [3]rotaxane plus cantilever is constrained to reside 
in the unstable region (II) of Fig.[6l with L = 70.0 A. Typical structures encountered in this regime 
(labels la-3a) are shown in Fig. [8l 
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FIG. 8: Structures observed in the unstable region of the puhing simulations. The numerical labels 
are employed in Fig. [7l 
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FIG. 9: Potential Uo{0 and entropic TS{^) contributions to the PMF [Eq. i^] of the [3]rotaxane 
along the end-to-end distance The results are averages of three different pulling simulations. 
The error bars correspond to twice the standard deviation obtained from a bootstrapping analysis. 
The solid lines result from a spline interpolation of the available data points. Typical structures 
(labels 1-7) are shown in Fig. O 
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FIG. 10: Probability density distribution of the radius of gyration p{R'^) for the [3]rotaxane when 
^ is fixed at 7.0 A (dots), 47 A (open circles) and 64.9 A (crosses). Note the bistability along Rg 
when is fixed at the concave region of the PMF. 
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FIG. 11: PMF as a function of the molecular extension ^ extracted from force measurements using 
WHAM for rotaxanes with different numbers of threaded rings. The error bars correspond to twice 
the standard deviation obtained from a bootstrapping analysis. 
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FIG. 12: Schematic variation of the Helmholtz potential of the molecule plus cantilever as a function 
of ^ for different extensions L suggested by the phenomenological observations around the region 
of mechanical instability. The labels I-III correspond to the different stability regions in Fig. [6] 
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FIG. 13: Left panel: dependence of the F-L isotherms on the cantilever spring constant k during 
the extension of the [3]rotaxane. Right panel: ratio between the thermal fluctuations and the 
average in the force measurements. The cantilever spring constants employed are all expressed in 
terms of the cantilever spring constant used in the pulling simulations ko = 1.1 pN/A. 
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FIG. 14: Dependence of the instability of the F-L isotherms on the cantilever spring constant k 
for the [3]rotaxane. The right panel shows the critical values of the force as a function of k. Here 
and F~ correspond to the values of the force when the F-L curves exhibit a maximum and 
minimum, respectively. The left panel shows the extension lengths that enclose the unstable 
region in the F-L isotherms, as well as the average end-to-end distance at the critical points 
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FIG. 15: Standard deviation in the force measurements ap = \J {F"^) l — ^ function of L 

for three different cantilever spring constants k during the pulhng of the [3]rotaxane. In each case, 
the dotted line indicates the value of ^Jk/ 13 which sets the hmit between the stable and unstable 
branches in the extension, see Eq. (j40p . 
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FIG. 16: Effective potential Ul{£,) = 4>{€) + for the molecule plus cantilever for different 

values of the extension L. In the panels, the open circles correspond to Ul{(,), the full circles to 
the PMF (^(^), and the solid lines to the cantilever potential Vl(^). The cantilever spring constant 
employed is the same as the one used in the pulling simulations presented in Sec. IIII A[ Note the 
bistability in the effective potential for L = 70 A. 
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FIG. 17: The upper panels show the probabiUty density distribution of the force measurements 
Pl{F) = {5[F — F{r,L)\)L (in pN~^) during the extension of the [3]rotaxane using cantilevers of 
varying stiffness. The force function F{r,L) is defined by Eq. ([3]). The lower panels show the 
associated spatial probability density distributions (in A^^) [Eq. (fTO|) ] for the [3]rotaxane 

plus cantilever. The color code is given in the far right. The spring constants are expressed in 
terms of /cq = 1.1 pN/A, the value employed in the pulling simulations presented in Sec. IIII A[ 
Note how bistability along ^, and hence blinks in the force measurements, arises for soft cantilevers 
and decays for stiff ones in accordance with Eq. (j42p 
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